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We apply a meron cluster algorithm to the XY spin chain, which describes 
a quantum rotor. This is a multi-cluster simulation supplemented by an im- 
proved estimator, which deals with objects of half-integer topological charge. 
This method is powerful enough to provide precise results for the model with 
a #-term — it is therefore one of the rare examples, where a system with 
a complex action can be solved numerically. In particular we measure the 
correlation length, as well as the topological and magnetic susceptibility. We 
discuss the algorithmic efficiency in view of the critical slowing down. Due to 
the excellent performance that we observe, it is strongly motivated to work 
on new applications of meron cluster algorithms in higher dimensions. 
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1 Introduction 



The functional integral formalism of quantum physics deals with infinite di- 
mensional integrals, which can only be computed explicitly in a few simple 
situations. A non-perturbative method to tackle models, which are not an- 
alytically soluble, starts by a regularisation to a finite number of degrees of 
freedom. This is usually achieved by a lattice discretisation of the time (in 
quantum mechanics) or of the space-time (in quantum field theory). In a 
finite volume the functional integral is then given by a finite set of single 
variable integrals. One tries to compute them numerically and — based on 
a variety of such results — to extrapolate to the continuum and to infinite 
volume. The transition to Euclidean time is very helpful to speed up the 
convergence of the integrals. 

However, the number of integrals still tends to be so large that straight 
numerical integration is hopeless. Instead one performs Monte Carlo simu- 
lations to generate a set of paths or field configurations with the Boltzmann 
probability distribution given by the Euclidean action. Thus one evaluates 
the expectation values of the observables of interest directly at finite inter- 
action strength (in contrast to perturbation theory). On the other hand, 
one has to face errors due to the limited statistics and uncertainties in the 
extrapolations. 

Hence it is essential to optimise the algorithmic tools for such simulations. 
The Metropolis algorithm is the most established procedure, but in many 
cases it is far from optimal. It can be refined to cluster algorithms [1, 2] 
which are by far more efficient for some set of models (for a review, see 
Ref. [3]). Unfortunately this set, where it could be applied successfully, is 
still quite small — in particular it excludes gauge theories up to nowJ§ But in 
the light of the striking success in specific spin models, it is highly motivated 
to explore cluster algorithms further. Here we present successful tests on 
(classical) spin chains, which describe quantum mechanical systemsJl 

Unlike the Metropolis algorithm, cluster algorithms do not proceed from 
one configuration to the next by updating single spins, but by flipping whole 

2 There have been proposals for cluster algorithms for U(l) gauge theory [4], but a 
breakthrough in the performance is still outstanding. For the treatment of a discrete 
gauge group, see e.g. Ref. [5]. 

3 The motivation we are giving here are functional integrals in quantum physics, but 
cluster algorithms have a much broader range of applicability, which also reaches out to 
fields like solid state physics and biology; for recent examples, see Refs. [5-8]. 
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clusters of them. First, this is promising in view of the thermalisation time 
needed in the beginning of a simulation. Later one expects to generate with 
the cluster algorithm well de-correlated configurations (which are needed for 
the measurements) with a modest number of simulation steps. In addition, 
multi-cluster algorithms can often be combined with an "improved estima- 
tor" , which allows for the inclusion of lots of configurations that do not need 
to be Monte Carlo generated explicitly. All these properties help to reduce 
the computer time required to measure an observable numerically to a given 
accuracy. This will be clearly confirmed for the system under consideration 
in this work, in particular as one approaches the continuum limit. 

We apply this technique to the 0(2) spin chain (the XY model), to be 
described in Section 2. There is a neat way to attach a half-integer topologi- 
cal charge to each cluster. This is the basis of a meron cluster simulation, 
which was first applied to a 2d 0(3) model on a triangular lattice (with a 
constrained maximal angle between neighbouring spins) [9], see also Ref. [10]. 
In this framework an improved estimator is extremely powerful. Variants of 
the meron cluster algorithm can also handle fermionic spin models success- 
fully [11]. 

In Sections 3 and 4 we present a novel application of this algorithm [12]. It 
enables us to approach the continuum limit much better than the Metropolis 
algorithm, and to suppress the notorious "critical slowing down". As in 
the original application, it is powerful enough to even explore the system 
with a 0-term. In almost all cases the simulation of a system with a complex 
Euclidean action is hardly feasible so far (perhaps up to a region of very small 
imaginary parts) Q Here we present one of the rare exceptions. Moreover, 
the constraint on the angles is not necessary in our case; the latter was 
required for technical reasons in the original application [9] (though it did 
not affect the universality class). Section 5 is dedicated to our conclusions 
and an outlook on potential applications of the techniques discussed here to 
a system of light quarks at high temperature. 

2 The quantum rotor 

We deal with a free scalar particle of mass m on a circle of radius 1, i.e. a 
quantum rotor. Its position is given by an angle <p(t), where t is the Euclidean 
time, so the Lagrangian reads L = ^Lpif) 2 . We consider the propagator G 

4 For reviews of the situation in QCD at finite baryon density, we refer to Refs. [13]. 
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between the end-points <p(0) = and <p(T) = (T > 0), i.e. we assume 
periodic boundary conditions. In the path integral formulation it can be 
decomposed into disjoint contributions with different winding numbers v. It 
is therefore the simplest quantum system with topological sectors. As in 
QCD we can also insert a #-term in this summation, which leads to 

+oo 

G(0,T;0,0) = Gv(0,T; 0,0) e ive 

-z/ 2 + ivO 



m 
2^f 



2mix 2 



ex p 



T 



v=—oo 

oo 



2tt ^ L 2m V 2tt 

v=— oo 

is the free propagator (27ri/, T|0, 0) on a line, and we set h = 1. It is 
sufficient to consider 9 G [0, 7r]. 

In QCD it appears natural that a #-term should occur, hence it is a mys- 
tery — known as the "strong CP problem" — why the observed #-angle is 
zero (or very close to it). In our case, a finite (9-angle does describe a physical 
situation, if we assume the particle to carry an electric charge q and a mag- 
netic flux $ to cross the circle. Then we identify 9 = 27rg$ (Aharonov-Bohm 
effect). Unlike QCD, we can evaluate in the present toy model the effect of 
9 > precisely, see Section 4. 

We now discretise the period T in L equal steps. We use lattice units, i.e. 
we set the step length T/L = 1. We denote y?(t = j) = ipj, and periodicity 
implies (fj+L = fj- The standard action on this temporal lattice reads 

L 

S s [(p] = m^J [1 - cos(^ j+ i - cpj)] . (2.2) 
j=i 

In contrast, the perfect lattice action [14], which is obtained from an infi- 
nite iteration of renormalisation group transformations^ distinguishes the 
topological sector that the particle may enter at time slice j 

L 



'For scalar particles in field theory this method is discussed in Refs. [15]. 
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In the functional integral all these sectors Vj e % are summed over, which 
reproduces the exact continuum result. 

This discrete system has another interpretation as a spin chain. On each 
site j a classical spin Sj = (S^, S^) of length Sf = 1 is attached. This is 
the XY model, which has a global 0(2) symmetry. If we stay with periodic 
boundary conditions and assume only nearest neighbour interactions with 
some coupling c, we arrive at the partition function 

L 

Z = Tr e-W] ? H ^s] = - c J2 • S, . (2.4) 

The trace means the sum over all spin configurations [S] = (Si, S 2 , ■ ■ ■ , £l), 
and (3 is an inverse temperature. If we identify the constants as (3c = m, we 
obtain the standard lattice path integral of the quantum rotor at 9 = (up 
to an additive constant in the action); the angle <pj describes the direction 

— * 

of the spin Sj, and also the spin model can be generalised by a #-term. 

The standard discretised system does not have natural topological sec- 
tors, because all configurations can be continuously deformed into one an- 
other (in contrast to the case of continuous time). Still one often introduces 
topologies, which is, however, ambiguous. The most obvious option is the 
geometric charge [16], which can be formulated analogously for instance in 
N- dimensional 0(N + 1) models, or in 4d Yang-Mills gauge theories [17]. In 
our case, the geometric charge amounts to 

1 L 

where Aipj = (<Pj + i — <Pj) G (— 7r, 7r]. 

We will also consider alternative formulations. 

2.1 Observables 

We are going to extract the correlation length £ as usual from the exponential 
decay of the connected 2-point function (resp. a cosh function due to the 
periodic boundary conditions). £ sets the scale of the system, and physically 
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sensible results usually require 

1<£«L. (2.6) 

The first (second) inequality implies that discretisation artifacts (finite size 
effects) are harmless. 

Our observables are the topological and the magnetic susceptibility, 

Xt = l((» 2 )-M 2 ) . ( 2 -7) 

1 L 
Xm = Z ((M 2 )-(M) 2 ) , M = J2Sj. (2.8) 



i=i 

Let us assume the cosh shape of the correlation functions to hold at all 
distances. Then Xm can be computed as follows, 



Xrn = l{(Y, 3 *y) = 1+ ib&fy 

3=1 3=2 

1 + ST [e-V-m + e -(x+i-i)/€j = 2 1 ~ _ J . (2 . 9) 



We now assume in addition the inequalities ( 12.61) to hold. In fact our simula- 
tions — to be presented below — were performed consistently^ at L/£ ~ 20. 
Thus the term e~ L ^ can be safely neglected, which leads to 



^2 + -^ + 0(r 3 ). (2.10) 



3 A meron cluster simulation of the XY model 
3.1 The algorithm 

We start by briefly reviewing the multi-cluster algorithm for O(N) models 
[2, 18], and in particular its extension to a meron cluster algorithm [9]. 

A step of the multi-cluster algorithm begins by building clusters, which 
are sets of neighbouring spins. To this end, a random direction r is chosen 



D We refer here to the correlation length £ at 8 = 0. 
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in an isotropic way (|r| = 1), and each spin Sj is split into Sj = (Sj ■ r) f and 

Sj- = Sj — Sj . A virtual bond is set between the sites j and j + 1 with the 
probability 

|l_ e -^ + , if ■ 3} : §} +1>0 (31) 
[ otherwise . 

Then a cluster is composed of neighbouring spins connected by bonds; 
it may also consist of a single spin, if the latter is disconnected. A step 
of the algorithm ends with "flipping" each cluster with the probability 1/2 
[1]. Flipping a cluster means that all its spins are mirrored at the plane 
perpendicular to r, Sj — > —Sj, Sj- — > Sj 1 . This algorithm respects 
ergodicity and detailed balance [2]. 

Let us numerate the clusters with c = 1, 2, . . . . A topological charge Q c 
can be assigned to each cluster based on the difference of the total topological 
charge v (i.e. the winding number) of the chain when the cluster is in its 
initial orientation, and after it has been flipped [9], 

q ^c— initial "c— nipped ^ ^ 

On the right-hand-side we use the geometric charge (12.51) . Q c is a half-integer, 
which remains unchanged if any other clusters are flipped, so it is determined 
locally [12]. To illustrate this important property, let Sj ... S}~ be the spins 
of a specific cluster. We denote the sum of the relative angles Aipj G (— tt, tt] 

(cf. eq. (12.51) ) between the successive neighbouring spins I to n as S[S n , and a 
flipped spin is written as S' t . The cluster charge only depends on its boundary 
spins and the two neighbouring spins of the adjacent clusters, 



Q Q I C C I C C Ql Ql Ql Q Q Ql 



(3.3) 



The charge is the same, regardless whether the neighbouring cluster is flipped 

or not, provided that Sj-±Sj — Sj-iS'j = S'j_ 1 Sj — Sj^S'j. In fact, this holds 
generally, which is easy to show by distinguishing different cases of the angles 

SjSj^i and SjS'j_ 1 . 

It is a peculiarity of the spin chain that there cannot be any loop inside 
a cluster, as we see from prescription (13. ip . Thus the cluster charges are 
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limited to the values 1/2, and —1/2, and the corresponding clusters are 
denoted as meron, neutral cluster and anti-meron, respectively. 

The property that Q c is determined locally for each cluster enables us to 
construct an improved estimator, which will be applied as a powerful tool in 
this work. With N c clusters, 2 Nc configurations can be obtained by cluster 
flips, which could enter the statistics (without the need for a Metropolis 
accept-reject step). In practice it is not optimal — or not even possible — to 
include all of them (we encountered N c values up to 40), unless this average 
can be evaluated analytically. Of course these configurations are not fully 
independent because they are all affiliated to the same direction r. 

3.2 Cluster statistics 

It is common lore that the "characteristic" cluster size follows the correlation 
length. Taking a close look at this property, we found that the statistical 
distribution of cluster with length s can be fitted well to a sum of three ex- 
ponentials, X^=i °i exp(— s/ Si). For large clusters the first term is dominant. 
Its decay is given by s\ = 1.000(1)£, in precise agreement with the expecta- 
tion. This supports the interpretation of the clusters as physical degrees of 
freedom, which is the basis of the meron picture employed here. At smaller 




20 40 60 80 100 

Cluster size 



Figure 1: The statistical distribution of the cluster sizes, measured at two dif- 
ferent correlation lengths. The density of large clusters decays exponentially, 
where the characteristic size coincides with the correlation length. 

cluster sizes the curve is steeper than the first exponential alone, see Figure [TJ 

7 In higher dimensions this vital property can only be achieved by imposing constraints 
on the maximal angles between neighbouring spins [9], cf. Section 1. 
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meron density 


cluster size s 


2.5 % 


0.401(5) £-0.10(3) 


12.5 % 


0.644(3) £ + 0.23(3) 


22.5 % 


1.16(1) £ + 0.93(7) 



Table 1: The cluster sizes corresponding to three specific meron densities. 

The leading sub-dominant exponential has a short range of S2 = 0.111 (2)£. 
We add that the continuum limit leads to a stable fraction of 58.0(1)% clus- 
ters of the minimal size s — 1. 

Next we consider the fraction of merons among the clusters. At large £ it 
amounts to 0.1581(5) £" a542{2) (of course the same holds for the anti-merons) . 
Obviously large clusters have a higher probability to carry topological charge. 
In the limit of a very large size s one finds 25% merons; around s = 2 £ one 
already arrives approximately at this asymptotic number. 

To provide an intuitive argument for this property, let us assume for 
instance a direction f— (0, 1), and we measure the spin angles (fj G (— 7r, it) 
relative to the x-axis. Again we consider some cluster with the spins Si . . . Sk 
and we assume <pi G (0, tt/2). Now the spin angles of this cluster describe a 
discrete path in (0,7r). For large £ (resp. large (3) the relative angles Aipj of 
adjacent spins are small. In particular we can assume the continuum limit 
to be approached to a point where the probability for any \Aipj\ > tt/2 is 
negligible. Hence for small clusters ifk is likely to be in (0,7r/2) as well, 
so that the cluster is neutral. However, in very large clusters (fi becomes 
irrelevant for the endpoint which can be in (0,7r/2) or in (tt/2,tt) with 
equal probability. This implies an equal number of neutral clusters and 
merons. Analogously, large clusters with ipi G (tt/2, it) are equally likely to 
be neutral or anti-merons. 

To quantify the increase of the meron density as s rises, we specified three 
densities and Table [1] displays the corresponding sizes s. 

3.3 Efficiency 

For comparison, we consider as a unit of computation time a process that 
could modify the whole chain: for the multi-cluster algorithm, this is what 
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we have described before as one algorithmic step; for Metropolis, it means 
one sweep to tackle each spin in the chain. We repeat that we performed 
our tests at a chain length L ~ 20 ■ £, so that finite size effects are strongly 
suppressed]^] 



-f 



Metropolis > 
fit : 7*exp(£/2.4) 
multi-cluster > 
fit:16T; 51 



Figure 2: The thermalisation time T t h er m for the energy as a function of the 
correlation length £. In the multi-cluster simulation r t h erm increases only 
slowly in £, which is in contrast to the Metropolis algorithm. The results are 
averaged over a variety of cold and hot starts. 

First we consider the thermalisation time with respect to the energy. 
For the Metropolis algorithm r t h erm grows exponentially with the correlation 
length, rtherm — 7(2) exp(£/2.4(2)) . On the other hand the data for the 
multi-cluster algorithm follow a power law, r t herm — 16(2)£ 0,51 ( 7 ) as shown in 
Figure [2l In particular this ensures a striking advantage at large £, when we 
approach the continuum limit. 

We proceed to the stage where the thermalisation is completed and we 
consider now the (exponential) auto-correlation time r a , again with respect 
to the energy. Hence the auto-correlation function is fitted with exp(— r/r a ), 
where r is the algorithmic time. The values of r a are plotted in Figure [3] (on 
the left) at various £ for the multi-cluster algorithm; we observe r a oc £ 7 with 
a dynamical critical exponent 7 = 0.52(3). 

For the Metropolis simulation, the data can be fitted well with a sum of 
two exponentials, cexp(— r/r i) + (1 — c) exp(— r/r a2 ), and Figure [3] (on the 
right) shows the corresponding results r a \ and t Q 2. The corresponding critical 

8 Efficiency studies in the two dimensional XY model, with multi-cluster and single 
cluster algorithms, were presented in Refs. [18,19]. 
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for Metropolis 
fit:1.11*E 1 - 



1 10 1 10 

£ E 

Figure 3: The auto- correlation time r a for the energy. The plot on the left 
shows the result obtained with the multi- cluster algorithm, where we observe 
a modest increases r a oc £ a52 . On the right we show the Metropolis result. It 
can be fitted with two exponentials, with a dominant exponent of 1.89. 

exponents amount to 1.89(6) and 0.89(4). What ultimately matters is the 
dominant exponent ~ 2, which is reminiscent of the random walk diffusion 
of local changes on the chain. 

For the (squared) topological charge (z/^) 2 (which is relevant for Xt), 
the growth of the auto-correlation time is exponential with Metropolis, as 
Figure H] shows. On the other hand, the auto-correlation practically vanishes 
with the multi-cluster algorithm. This de-correlation is due to the large- 
scale changes performed on the chain. The multi-cluster algorithm reveals 
here most clearly its potential in overcoming the critical slowing down. 

4 Results for the observables 
4.1 Correlation length £ 

In Subsection 3.2 the correlation length (at 8 = 0) has been anticipated. 
We now consider its relation to the inverse temperature (3. The numerical 
and theoretical [14] results with the standard action match perfectly, as the 
plot in Figure [5] on the left shows. On the right we add results for £ at 
non-zero vacuum angles 8. For increasing j3 the correlation length of the 
standard action approaches the perfect action value, £/(2/3) = 1/(1 — 9/ir). 
A divergence at 8 —> 7r has also been observed in the 2d 0(3) model [9]. 
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Metropolis data ^ 
fit: 0.51 *exp©1.31) 



Figure 4: The auto- correlation time r a with respect to the topological charge 
squared. The plot shows its exponential increase with £ for the Metropolis al- 
gorithm. In contrast, it practically vanishes with the multi-cluster algorithm. 
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Figure 5: The correlation length as a function of (3 for the standard action. 
The simulation results at 9 = (on the left) are in excellent agreement with 
the theoretical prediction [14]- At finite 9 (shown on the right) and increasing 
f3 they converge towards the perfect action result, which describes the system 
in the continuum. 



4.2 Topological susceptibility xt 

Figure [6] shows the accurate agreement of the measured topological suscep- 
tibility (eq. (12. 7p ) at 9 = with the theoretical formula of Ref. [14]. 

To calculate Xt{@) a l so at 6* 7^ 0, the probability p(u) is needed to an 
extremely high precision. It proves to be Gaussian, see Figure [7] (on the 
left). An improved estimator can be used here which captures all cluster 
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Figure 6: The numerical results for the topological susceptibility Xt (at 9 = 0), 
in precise agreement with the theoretical prediction for the standard lattice 
action [14]. The result for the perfect action coincides with the continuum 
susceptibility at any correlation length. 




orientations by simple combinatorics. This enhances the statistics drastically, 
and it allows us to reach probabilities of 0(1CT 12 ) with only one million really 
generated configurations. Subsequently we can use an analytic expression for 
the distribution p(u). 



1 

0.01 
e-04 
e-06 
e-08 
e-10 
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simulation : p = 1, L = 60, 1e6 configurations 
0.2801 (2)'exp(-1 4.784(8) g (v "2) / L ) 




Figure 7: On the left: the probability distribution for the topological charges, 
p{v), which follows a Gaussian. An improved estimator is mandatory to cap- 
ture even tiny probabilities. On the right: Xt(6)£ a t different vacuum angles, 
against the length L of the spin chain. At large L it becomes independent of 
9 ( except for 9 = ir). 

The results for the ^-dependent susceptibility Xt(6) = \ (( z/2 )e ~ ( u )e) 



13 



were obtained by relying on the Gaussian distribution p(u) that we identified. 
Note that both terms contribute. Xtifi) is rea l due to the parity symmetry, 
which implies the symmetry in the sign of v. Figure [7] (on the right) shows 
its dependence on L for various values of 9. (To set the scale, we still refer 
to £ at 6 = 0.) It converges for large L to the value of Xt(0), for any 6 ^ n. 
This convergence slows down as 6 rises, and it collapses at 9 = it. 

Such precise results for a system with a complex action are very rare; other 
examples with a #-term were obtained with the meron cluster algorithm for 
the 2d 0(3) model [9] and for an SU(N) quantum spin ladder [20] 

4.3 The topological susceptibility from cooling 

For comparison we also consider this susceptibility based on topological 
charges obtained from "cooling" [22]; we denote it as Xt.cooi- To this end, 
the chain is smoothed before a measurement: a spin is chosen at random 
and rotated so that the action is minimised. This process is iterated until it 
converges. In lattice gauge theory a long cooling process (for the plaquettes) 
ultimately leads to the trivial configuration, so one tries to read off a topolog- 
ical charge from some earlier plateau (for instance by monitoring the energy). 
Here the situation is simpler because the cooled configuration stabilises, so 
this sense) unambiguous. The question remains how it is related 
to the original configuration, which has the correct statistical weight (unlike 
the cooled configuration). We may compare the cooled charge to the original 
geometrical charge u^ g \ with the obvious inequality |z/( cool )| < |i/ g )|. The 
ultimate criterion is how well Xt.cooi approximates the continuum value of %t- 
The result is plotted Figure The convergence of Xt.cooi to the continuum 
limit is just as quick as for the geometrical charge without cooling. 

4.4 Magnetic susceptibility Xm 

Figured (on the left) shows our results for the magnetic susceptibility Xm, see 
eq. (12.81) . They are in excellent agreement with the approximation given in 
eq. fl2.10p . We mention that Xm corresponds to the mean cluster size when 
one employs the single cluster algorithm [2], which favours larger clusters 
than the multi-cluster method. 

9 For alternative methods to simulate models with a 6>-term, see e.g. Refs. [21]. 
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Figure 8: The topological susceptibility obtained from "cooling". For increas- 
ing £ converges to the continuum value with the same speed as the use of 
the geometrical charge on the uncooled configurations. 



2.16 
2.14 
2.12 
2.1 
2.08 
2.06 
2.04 
2.02 
2 

1.98 



continuum 
simulation 
2 + 1/(6 £?) 



10 



15 



20 




Figure 9: The numerical results for the magnetic susceptibility. On the left 
we show the result at 9 = 0, which follows the formula A2.10)) to an excellent 
accuracy. The plot on the right is our result for Xm{9)/£,(9). The data are 
reliable over the plotted range, before getting too close to tt. This result was 
accessible with the multi-cluster algorithm along with an improved estimator. 



At non-zero 9, the result emerges from re-weighting the data of 9 = 
according to the prescription in Ref. [23]. The results are shown in Figure 
[9] on the right. We applied an improved estimator by flipping a number of 
clusters. Here we do consider the ^-dependence of the correlation length, 
which is used to build the dimensionless ratio plotted in Figure It grows 
rapidly as 9 approaches tt, cf. Figure hence we involved quite long spin 
chains. Once £(9) exceeds L, the considered ratio drops towards 0, as we see 
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from eq. (j2.9p . However, for a continuum limit 9 — > 7r which respects L 3> £ 
the plateau at relatively large 6* suggests a value Xm{&) ~ 1- 

In one case, (3 = 5, L = 300, we indicate jackknife errors. In the other 
cases they are similar, i.e. again very small up to the close vicinity of 9 = it. 
This is very remarkable in view of the generic difficulties to obtain numerical 
results for models with a significant imaginary part in the Euclidean action. 

The statistics includes several millions of configurations, and it is still 
enhanced thanks to the improved estimator. In particular, statistics in a 
fixed topological sector can be cumulated by flipping neutral clusters. In 
this way it could be simulated within a few weeks on a 2 GHz machine; with 
the Metropolis algorithm this measurement is hardly feasible. 

5 Conclusions and outlook 

We presented a new and successful application of the meron cluster algorithm. 
In the spin chains that we considered, it provides precise simulation results 
in a highly efficient way. We assigned a half-integer topological charge to 
each cluster, which is the basis of a powerful improved estimator. This yields 
accurate result for the present toy model with a 0-term — an issue, which is 
still outstanding in QCD. 

Here it was possible to suppress the notorious problem of critical slow- 
ing down. With respect to the topological charge we overcome this problem 
completely. Regarding the energy, the dynamical critical exponent is reduced 
by almost a factor of 4. 

This provides a strong motivation to search for applications of this tech- 
nique also in higher dimensions, i.e. in field theoretic models. In particular 
an application of the meron cluster algorithm in the 3d 0(4) model appears 
promising — that model can be interpreted as an effective description of 
QCD with two light quark flavours at high temperature. 

We add a very rough estimate about the feasibility of that project. If 
£ (9) behaves similarly to Figure we control the finite size effects quite well 
up to 9 fa 0.9 tt for instance with L = 32. Compared to the spin chains that 
we considered to measure Xm(9)/£,(9), this implies a factor of 0(100) for the 
number of lattice sites (along with a factor of 6 for the generators of the 
symmetry group). Tiny error bars as we obtained in Figure [9] (on the right) 
may be relaxed without problems, say by a factor >3, so that the required 
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statistics decreases by an order of magnitude. Comparing now to the compu- 
tational effort which was necessary in d = 1 (cf. last paragraph in Section 4), 
and considering the option to use a number of processors simultaneously, we 
estimate that the 3d 0(4) model at finite 9 can be solved to a good precision 
with the meron cluster algorithm within less than one year. 
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